Galactic halos of self-interacting dark matter 
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Recent, very accurate simulations of galaxy formation have revealed that the standard cold dark 
matter model has great difficulty in explaining the detailed structure of galaxies. One of the major 
problems is that galactic halos are too centrally concentrated. Dark matter self-interactions have 
been proposed as a possible means of resolving this inconsistency. Here, we investigate quantitatively 
the effect of dark matter self interactions on formation of galactic halos. Our numerical framework is 
extremely simple, while still keeping the essential physics. We confirm that strongly self-interacting 
dark matter leads to less centrally concentrated structures. Interestingly, we find that for a range of 
different interaction strengths, the dark matter halos are unstable to particle ejection on a timescale 
comparable to the Hubble time. 

PACS numbers: 95.35. +d, 98.62. Gq, 14.80.-j 



I. INTRODUCTION 

The concept of dark matter was originally introduced 
by Zwicky in 1933 |Q. to explain the behaviour of indi- 
vidual galaxy clusters. Since then the dark matter has 
been determined to be essential for explaining a vast 
number of astronomical observations. The cold dark 
matter model has, with modifications, been very success- 
ful in explaining how structure forms in the universe [Q. 
However, it was recently realised that the CDM model is 
apparently unable to explain the dynamics of individual 
galaxies. Very high resolution simulations have shown 
that there should be much more sub-structure in the ha- 
los of typical galaxies than is observed. Numerical sim- 
ulations predict that the local group should at present 
contain about 1000 distinct dark matter halos. Observa- 
tions yield a number which is a factor of ten lower 0,0. 
Also, lensing measurements of clusters seem to indicate 
a constant density core, in strong contrast with numeri- 
cal simulations ||||. Also, the central density profile in 
galactic halos is predicted to be much steeper than in real 
galaxies |o]-|lC|. However, this is not a completely settled 
issue. For instance Kravtsov et al. II find that CDM 
simulations are consistent with observations. 

If the problem with CDM halos persists, then it could 
possibly be remedied if structure formation is somehow 
suppressed on small scales. One such possibility is that 
perhaps the dark matter is not cold, but rather "warm" , 
i.e. dark matter particles have masses around 1 keV and 
therefore have significant thermal motion around the 
time of matter-radiation decoupling |T^ ]. That would 
suppress structure formation on galactic scales and be- 
low. However, warm dark matter may have problems in 
describing properly the properties of clusters. Another 
possibility is that the initial power spectrum has a cut-off 
at some wavelength, corresponding to the substructures 
in galaxies , so that no structure below this scale will 
grow initially. 

A quite different possibility is that the discrepancy has 



to do with the interaction properties of the dark mat- 
ter, not just the mass of the individual particles. It 
was recently proposed by Spergel and Steinhardt fllij , 
that dark matter with strong self-interactions could make 
the cold dark matter model consistent with observations. 
Strongly self-interacting dark matter will tend to pro- 
duce halos with shallower core density profiles, thereby 
alleviating the problem of the central mass concentration. 
Although the simple arguments provided in Ref. |l4| are 
quite convincing, it seems very important to test the im- 
plications of self interacting dark matter on a more quan- 
titative basis. That is the purpose of the present paper. 
With the help of a very simple numerical scheme, we 
solve the Boltzmann equation describing the phase-space 
evolution of collisional dark matter and are thus able to 
calculate final state density and velocity distributions. 

Intriguingly, we find that there is a range of param- 
eters for which there exist no stable equilibria, even on 
a relatively short timescale. In general, however, we are 
able to confirm the predictions of Ref. |l4|] , namely that 
strongly self interacting dark matter produces halos with 
shallower density profiles and thus can remedy the prob- 
lems that the standard cold dark matter model has in 
explaining galactic dynamics. It should be noted here 
that self interacting dark matter has been considered pre- 
viously, in order to study possible effects on the initial 
linear power spectrum Jl5]-|l7j. However, the type of self 
interactions investigated in the present paper will have 
no discernible effects on the initial power spectrum. 



II. THE PHYSICS OF SELF-INTERACTING 
DARK MATTER 



In the standard cold dark matter model, the dark mat- 
ter particles were once in thermal equilibrium in the early 
universe. As the temperature dropped below their rest 
mass, however, their abundance was exponentially sup- 
pressed. Eventually annihilation reactions were no longer 
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effective and the abundance of dark matter particles froze 
out at some particular value [T^] . In this picture, the 
dark matter particles are presently extremely weakly in- 
teracting and halos consisting of them behave as colli- 
sionless systems. However, dark matter particles could 
in principle have very strong self-interaction, as long as 
their reaction cross section with standard model parti- 
cles is small. That is for instance the case with the non- 
thermally produced particles proposed by Chung, Kolb 
and Riotto Jl8| ] . In the present paper we shall assume 
only that dark matter particles are completely decoupled 
from all other species. The dark matter self-interactions 
are modelled as point-like, i.e. the interactions are medi- 
ated by some very massive boson. In that case, the cross 
section for two-body scattering, cr, is energy-independent. 
This means that the reaction rate for a dark matter par- 
ticle moving in a dark matter halo is 



d 3 p 



fpVrelCT, 



(1) 



where u re i is the relative velocity between the incident 
particle and the target and f p is the phase-space dis- 
tribution of dark matter particles. As will be explained 
later, the cross-section that divides the very weakly inter- 
acting regime from the strongly interacting one is given 
roughly by 
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which of course is many orders of magnitude larger than 
the interaction cross section for standard cold dark mat- 
ter. 



III. NUMERICAL SCHEME 

The fundamental equation describing the phase-space 
evolution of self-interacting dark matter is the Boltzmann 
equation. We shall in the present paper neglect the Hub- 
ble expansion of the universe. In that case, the Boltz- 
mann equation takes on the form 



dt dt V ' dr 



ov 



(3) 



where the right-hand side is a collision operator describ- 
ing possible scattering/annihilation reactions. 

For normal, collisionless matter, the right hand side is 
zero, but in our case it is a phase space integral of the 
reaction matrix element. 

In order to keep the numerical work quite simple, while 
still retrieving the essential physics, we shall assume 
spherical symmetry of /. In that case, there are only 
four independent phase space coordinates, (t, r, v r , v±)- 
Instead of v± we use the specific angular momentum, 
j = v±r, which is a conserved quantity for collisionless 
systems [^0| . The left hand side of the Boltzmann equa- 
tion is then |E0(1 
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GM (r) 



f- (4) 



Our numerical scheme for solving the Boltzmann equa- 
tion is somewhat similar to that of Rasio, Shapiro and 
Teukolsky H]. We mimic the phase space distribution 
by particles, distributed according to the initial /. Then 
these particles are moved in phase-space according to the 
Boltzmann equation. Thus, the force on each phase- 
space "particle" is calculated using the mean gravita- 
tional field produced by all particles, GM (r) jr. All mean 
field quantities, like M(r), are calculated on an Eulerian 
grid in phase-space and interpolated using spline interpo- 
lation. Thus, the code is a hybrid between a real Boltz- 
mann code and an N-body code. This is a standard way 
of solving the collisionless Boltzmann equation, both in 
gravitational dynamics and in plasma physics p2| , p3| . 

The right-hand side depends on the specific scattering 
process. In the case of two-body scattering, it can be 
written as |l8| 

C[f] = ^J d 3 p 2 d 3 p 3 d 3 p 4 A(f, f 2 , f 3 , U) (5) 

X E Nl2^34^(Pl + P2 - P3 ~ P4)(2^) 4 . 

Here we have used the definition d 3 p = d 3 p/ {(2tt) 3 2 E). 

\M\ 2 is the interaction matrix element and pi is the 
four- momentum of particle i. The phase-space factor is 
defined as A = / 3 / 4 (l - /)(1 - h) - // 2 (1 - /a)(l - h). 

However, since we are considering phase space to be 
sampled by a distribution of particles, changing the phase 
space density due to collisions amounts to letting the 
particles scatter on each other. 

The interaction rate for a given particle can therefore 
be written as in Eq. (Q) 



r = 



d 3 p 
(2^)3 



fpVrelCT ~ ncr^disp, 



(0) 



where fdisp is the velocity dispersion of the system. 
Therefore, in this simple approximation, the interaction 
rate of a dark matter particle is given in terms of the 
simple mean field quantities n and v<xi S p- 

The above system of equations can be relatively easily 
solved for given initial conditions, and quantities such 
as the radial density distribution and velocity dispersion 
can be recovered. 



IV. RESULTS 

We have solved the system of equations described in 
the previous section for different values of the interac- 
tion cross section, a. The Boltzmann equation is non- 
dimensionalised in the following way: Given a total mass 
Mo of the system and a characteristic radius Rq , the sys- 
tem is non-dimensional ised by s olving in terms of the 
dynamical time r = \/Rf)/ GMq and the characteristic 
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velocity vq — y/ GMq / Ro . The interaction cross section 
should then be cast in units of ao = mR 2 /Mo, where m 
is the mass of the dark matter particle. Note that for 
o~ = 00, the mean free path in a typical system should be 
A ~ Ro, so that that ao is a natural dividing line between 
slow and fast interactions. 

In all the simulations we start with a homogeneous 
sphere of radius Ro- The velocity dispersion is assumed 
to be Maxwellian so that 

fcxexp(-v 2 /2v 2 -j 2 /2r 2 v 2 ). (7) 

The velocity dispersion, Vd, is chosen to be Vd = O.l^o- 
If all particles start at rest, they will pass through the 
center at the same time, and the simulation becomes un- 
stable |24|. However, letting the particles have a small 
initial velocity does not alter the final state of the whole 
system very much. It only results in the collisionless sys- 
tems being slightly less centrally concentrated. 

The system is allowed to evolve until t = 30r, in order 
to test stability of the code. 

A. Collisionless systems 

The first simulation was done for a normal collisionless 
dark matter halo. Fig. 1 shows how the mass distribu- 
tion M(r) evolves with time. As expected, the system 
very quickly approaches equilibrium, essentially in just 
one dynamical time. The core profile of the relaxed sys- 
tem approximately follows a power law, p oc r~ 4 / 3 . Note 
that if the initial velocity of all particles was zero, the 
final system would have a density profile with p oc r -9 / 8 
instead p3]. Thus, our assumption of a small initial ve- 
locity leads to a less centrally concentrated halo struc- 
ture. However, since our purpose is to compare halos of 
dark matter with different interaction strengths, and not 
so much to make exact calculations, this assumption docs 
not matter much. 
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FIG. 1. The mass distribution, M(r), for the collisionless 
dark matter collapse at three different times. The solid line 
is for t — t, the dashed for t = lOr, and the dot-dashed for 
t = 20t. 



B. Collisional systems 

Ejection — As described above, a = ao divides the 
weakly and strongly interacting regimes. For a <C ctq, 
one would expect an evolution very similar to the fully 
collisionless system, whereas for a 3> a , the system will 
behave like a collisional gas, i.e. shocks can be produced 
and single particles will perform random walks in phase 
space. The halo should therefore behave like an entirely 
hydrodynamical core surrounded by a collisionless sys- 
tem. However, at a ~ ao the evolution is less obvious. 
As also noted in Ref. |Q, most self gravitating systems 
can eject particles. In normal collisionless systems, this 
can happen for two reasons, either a particle experiences 
a single close encounter which leaves it with positive en- 
ergy Hj| or it suffers many weak collisions, gradually 
increasing its energy |2(|. The first process is usually 
referred to as ejection, whereas the second is called evap- 
oration. This process makes the entire system unstable 
in the long term, since the only stable solution is a single 
pair of particles in a Kepler orbit, with all other particles 
being at infinity p7[ . However, for a normal galaxy of 
stellar objects, the instability time is vastly larger than 
the Hubble time, and for a collisionless dark matter halo 
it is for all practical purposes infinite. 

However, for a collisional halo, this need not be the 
case. Since particle interactions are point-like, a given 
particle can normally only gain positive energy because 
of a single scattering event, not because of many weak 
collisions. We can estimate the instability time in the 
following way: In a homogeneous mass distribution with 
a typical velocity vo = y GMq / Ro , a particle passing 
through the entire system has the scattering probability 

P ~ 1 - e - Ro/x , (8) 

where 



If the particle scatters, and acquires a positive energy, 
then it should not scatter again on its way out, if it is 
to be ejected. Therefore the total probability for the 
particle to be ejected is roughly 

Ejection Si (1 - e- Ro ' x )e- R °/ x P{E^ aX > 0). (10) 

-P(-Efmai > 0) is the probability that the energy after 
scattering is larger than zero. Normally this probability 
will be .P(-Efinai > 0) ~ 0.1 — 0.2, so that it is not negligi- 
ble. The probability is maximal when A ~ i? , as could 
be expected. The typical timescale for instability is then 
roughly 
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^instability — T J ^"ejection? 



(ii) 



From these extremely crude estimates one finds that 
the minimal possible instability time of of the order 
^instab iiity — few x r. Even if interactions are very strong, 
the system is unstable in the long term, because particle 
will escape from the system via diffusion |14| . However, 
the timescale for this diffusion process is much longer 
than the Hubble time. 
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FIG. 2. The mass distribution as a function of r at three 
different times and three different values of a. The solid lines 
are for t — r, the long-dashed for t — 15r, the dashed for 
t = 20r and the dotted for t = 25r. 

Numerical simulations — We have performed simu- 
lations for different values of er/c% specifically a/cro = 
1CT 3 , 1CT 2 , 1CT 1 , 1, 10, fOO, lOOO. Fig. 2 shows the evo- 
lution of M(r) for <j/<t = IfT 3 , I and 10 3 . For 
<t/o-q = I0 -3 , it is practically indistinguishable from the 
collisionless would be expected. For the very 

strongly interacting case a/o~Q = 1000, we find a very 
different behaviour. The system again settles into equi- 
librium very fast, roughly at t ~ few x r, but we also see 
that the core is much less dense because low entropy ma- 
terial is ejected. This is the result predicted by Ref. j jTjj . 

The intermediate case, a/ao = 1, never settles into a 
true long-term equilibrium. It starts out by approaching 
the collisionless equilibrium, but scattering interaction 



drives the mass distribution towards the equilibrium dis- 
tribution for the strongly interacting particles. However, 
on the same timescale the system loses particles due to 
the ejection mechanism described above, so that no true 
equilibrium is ever reached. The timescale for particle 
loss indeed seems to be ~ 10 x r. 

Observations — We cannot expect our simulations to 
fit observational rotation curves because of the approx- 
imations we have made, i.e. Hubble expansion was ne- 
glected and spherical symmetry was assumed. Neverthe- 
less, it is interesting to compare our results to observa- 
tions. In Fig. 3 we show the predicted rotation curves at 
t = 15r for the three different values, a = 10 -3 , I and 
I0 3 . The rotation curve for the strongly interacting halo 
is much less centrally peaked than the collisionless one. 

The crosses in Fig. 3 show the measured rotation curve 
for the typical low surface brightness galaxy UGCI28 ||. 
The observational rotation curve has an asymptotic cir- 
cular velocity of v = 200 km/s at R = 100 kpc. Tak- 
ing iio = 200 km/s and Ro — 100 kpc gives a mass 
of Mo = 1.35 x 10 12 Mq and a dynamical timescale of 
t = 4.7 x 10 s y. These numbers are typical for such 
galaxies. Notice that for this choice of mass and radius, 
oo = 8.5 x 10~ 23 cm, 2 mQ eV , i.e. corresponding to the <7o 
quoted in Eq. (ph. 
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FIG. 3. The rotational velocity at t = 20r for three dif- 
ferent values of a. The crosses are the observational rotation 
curve for the low surface brightness galaxy UGC128 M. 
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The rotation curve for the strongly interacting system 
provides a better fit to the observational rotation curve 
than that for the collisionless system. However, in all 
cases, the fits are quite poor. As mentioned, this should 
not be taken too seriously since our model involves some 
essential approximations. 



V. DISCUSSION 

We have performed quantitative calculations of how 
dark matter halos form in models with self interacting 
dark matter. For simplicity we assumed spherical sym- 
metry, neglected the Hubble expansion, and used a very 
simple prescription for the dark matter self-interaction. 

Our results have essentially confirmed the estimates of 
Ref. jL4|, in that dark matter halos with sufficient self 
interaction will have much shallower density profiles that 
normal collisionless halos. Thus, they can provide better 
fits to observational galactic rotation curves, while still 
being consistent with all other known data. 

A very curious feature is that halos with intermedi- 
ate self interaction (A ~ Rq) are unstable on a rather 
short timescale because they eject particles continuously. 
This effect was estimated to be unimportant in Ref. [jl4j , 
but our numerical simulations indicate that there exist 
regions of parameter space where the instability would 
probably have shown up in present day halos. 

As mentioned, our calculations assume spherical sym- 
metry and neglect the Hubble expansion. It will be very 
interesting to investigate the effect of self interacting dark 
matter using detailed N-body calculations ]2Sj] . Most 
likely, such simulations will confirm th e g eneral state- 
ments made by Spergel and Steinhardt nM as well as in 
the present paper. 
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